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ABSTRACT 


The spline fit curve is a convenient method for fitting a curve through a given set of 
points. This report describes a FORTRAN IV computer program which will calculate 
the spline fit curve, with function values, first and second derivatives, and curvature at 
any desired interpolated 1 points. 


FORTRAN PROGRAM FOR SPLINE FIT CURVE 
by Theodore Katsan Is 
Lewis Research Center 

SUMMARY 

The spline fit curve is a convenient method for fitting a curve through a given set of 
points. This report describes a FORTRAN IV computer program which will calculate 
the spline fit curve, with function values, first and second derivatives, and curvature at 
any desired interpolated points. , , ' ' 
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• .INTRODUCTION , . ■ 

If a set of function values corresponding to a set of arguments is given, there are 
several ways a curve can be fitted through these values so as to approximate the original 
function with these values. The classical way is by an n in degree polynomial for 
n + 1 points. However, this may not be satisfactory for a large number of points. 

Least squares approximations are satisfactory for smoothing, but if the function values 
are accurate, it is not desired to change these values as a least squares fit would. 

Another technique is to use fewer points as in four-point Lagrangian interpolation. But 
this does not lead to a smooth curve since the derivatives will not be continuous. 

A method that has received much attention«lately is the piecewise cubic, with con- 
tinuous first and second derivatives, commonly referred to as a spline fit curve (ref. 1). 
The spline fit curve is a mathematical expression for the shape taken by an idealized 
spline (thin wood or metal strip) passing through the given points. The spline fit curve 
gives a simple method of determining an approximating analytical curve which can be used 
in place of the original curve for interpolation, determining first and second derivatives, 
curvature, or integration. 

The spline fit curve is ideally suited for computer vise and has been successfully 
used in a number of computer programs. This report gives a FORTRAN IV program 
for the purpose of calculating the spline curve and its integral, with interpolated values 
of the function, first and second derivatives, and curvatures. This program has READ 



and WRITE statements for use as an independent program. The program can also be 
used as a subroutine with the input and output variables as subroutine arguments. 


METHOD AND THEORY 

Suppose that an interval a s x < b is divided into N intervals by a = x Q < < . . . 

< x N = b. The function F(x) is given at these N+l values of x by F k = F(x k ). Join 
each given point by a cubic polynomial and require that the first and second derivatives 
be continuous at these points. Each cubic has four coefficients for a total of 4N unknown 
coefficients. Each cubic is required to pass through the given points providing 2N con- 
ditions to be satisfied. The requirement of continuous first and second derivatives at 
the N - 1 interior points provides another 2N - 2 conditions, for a total of 4N - 2 con- 
ditions to be satisfied. Two additional conditions are necessary to determine the cubic 
spline curve uniquely. These are usually given as arbitrary end conditions. The end 
condition used in this program is that the second derivative at an end point is one-half 
the second derivative at the next point. 

The theoretical reason for using a cubic spline curve is that the cubic spline has the 
minimum value for 



of all curves passing through the spline points. This is proven in ref. 1. Since 



dx 


is nearly proportional to the strain energy of a thin, uniform spline with a small slope, 
the cubic spline fit will approximate the shape taken by a thin physical spline passing 
through the given points. The end condition of making the second derivative at the end 
point one -half the second derivative at the second point is equivalent to bending the spline 
beyond the last point slightly, instead of just allowing it to be straight. 

In reference 1 there is an extensive discussion of the mathematical theory of spline 
fit curves. Equation (4) of reference 1 is the basic equation to be solved. A complete 
set of equations is obtained by specifying values for the parameters X and X N of 
reference 1 (following eq. (4)). The values of X Q and X N are determined by the end 
conditions. The end conditions mentioned above result in X Q = X N = 1/2. The resulting 
complete set of equations can be expressed as a tridiagonal matrix equation which is 
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strongly diagonally dominant. The equation is solved by Gauss elimination. The equa- 
tion of the spline curve is then equation (3) of reference 1, with the derivative given by 
equation (2) and the second derivative by equation (1) of reference 1. The curvature K 
is calculated by 


K = 


[l + (F’) 2 ] 


372 


The radius of curvature is the reciprocal of the curvature. 


USE OF PROGRAM 

The spline fit curve can give very accurate values for first and second derivatives. 
To do this, however, the spline points must be accurately given. The derivatives are 
extremely sensitive to small errors in coordinates when points are close together. For 
this reason it is best to use as few points as is reasonable. If the curvature is slight, 
the spline points can be widely spaced. The greater the curvature, the closer the spline 
points should be (analogous to a physical spline). With a very sharp bend, the points 
should be quite close together and must be very accurate. 

If the curve has a large slope, the spline approximation is not as good. There is 
no trouble if the slope is less than one in absolute value. Much higher slopes can be 
handled satisfactorily, but more points are required. Great care must be used when 
using curves approaching the vertical within a few degrees, especially if there is appre- 
ciable curvature. 


Curves Given by Parametric Equations 

The spline curve fit is limited to curves of the form y = f(x). Thus, a circle could 
not be approximated by a spline fit curve. However, a curve of any shape can be ex- 
pressed parametrically (ref. 2). For example, x and y can both be given as a function 
of arc length s. That is, 


x = x(s) 
y = y(s) 

For example, for a circle of radius r with center at the origin, we have 
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X = 

y = 


r cos 
r sin 



Then the derivative dy/dx can be calculated by 


dy 

dy _ ds 
dx dx 
ds 

2 2 

The second derivative d y/dx can be calculated by 


dx 2 


dx 

ds 


a 

ds 2 


dy 

ds 


A 

ds 2 




The calculation of the first and second derivatives for parametric curves has not been 
included in the program. 


Use as a Subroutine 

This program can be used as a subroutine. Most likely, the READ and WRITE 
statements would be eliminated. The spline input variables (X, Y, and N) would be 
input arguments as well as the interpolation data (Z and MAX) if interpolation is 
done. The desired output variables would be the output arguments. Of course, 
calculation of unnecessary variables could be eliminated. 


DESCRIPTION OF INPUT AND OUTPUT 

Input. - The input data sheet is shown in figure 1. The quantities filled in are 
for a sample case. The output for this case is presented in the next section. N and 
MAX are both integers (no decimal point) and must be right -adjusted. The remaining 
input is real numbers (punch decimal point) in a 10-column field. 
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Figure L - Input form. 


The input variables are as follows: 


N 

MAX 


X - COORDINATES 
(SPLINE POINTS) 

Y - COORDINATES 
(SPLINES POINTS) 

X - COORDINATES 

(INTERPOLATION POINTS) 

Output. - An example of the 


Number of spline points (minimum 2, maximum 50) 

Number of interpolation points (minimum 0, maxi- 
mum 50) 

Set of X coordinates of the N spline points 

Set of Y coordinates of the N spline points 

Set of X coordinates where output data will be 
given 

output is given in figure 2. The first part lists the 


SINE CURVE 


LIST IE SPLINE 
X 

0 

0.52369680 
1.04719760 
1.67C79o30 
2. 09439510 
. 2.61799 389 

3. 14159271 


IIMPUT) IN 
Y 

0 

0. 50000000 
0. 86602540 
1. 00000000 
0. 86602540 
0. 50000UUO 
0 


CODROlNAIcS 


YALJE OF THE 
INTEGRAL OF 1 *D< 
3 

3.13494445 
3. 53:366853 
1.33356283 
l. 53365725 
1 . 85o38 L 33 
2.331 32576 


CUORl) HATES • 0FUVATIVES A NO CJKVATURtS AT SELECTED POINTS (NO. 


<< INPUT ) 

0 

0.50000000 
l.OOOOOUOU 
1.57079630 
7.00000000 
3. 00 COO 000 
3. 141 69271 
0. 52369880 


Y 1 1 NTERPJLA TE D1 
0 

0.4 7956131 
0. 8413386b 
1.00000000 
0. 90936258 
0.14336904 
0 

0. 50000000 


DY/OX 
L. 0336151 7 
0.8673919? 
3.54372633 
0. 6307330 iE -07 
-3 . 41 58 395 7 
-3.99733205 
-1.03361511 
0.85 66 7263 


OF PDHTS = 
027 7 0X2 
-0. 22541 757 
-3.44367548 
-3.86253815 
-1.01437925 
-0.92334693 
-3.29637536 
-0.22541 749 
-0.45393514 


:j*VA Tlrtr 
0. 75779933E-01 
0. 19005273 
0. 58486097 
1.0143 7925 
0.72637697 
0. 10164735 
3. 75778910E-01 
0. 19749248 


R A 31 US 0= Z UR V A T JR 
-13.1952831 
-5. 261 59755 
-1. 70913334 
-0. 93532453 
-1.373T435? 

-9. 83793373 
-13.1962837 
-5. 35343377 


Figured - Sample output. 
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spline points given as input, together with 



Y dX 


for I = 1, 2, . . . , N. 

The next output lists the X coordinates given as interpolation points in the input, 
with the interpolated values of Y, the first and second derivatives, the curvature, and 
the radius of curvature. 

The message OUT OF RANGE X is printed if there is extrapolation of more 

than one-tenth of either of the end intervals. Extrapolated values must be used with 
caution. 


COMPUTER PROGRAM 


Program Listing 

i) IMFMS ION X l 50 ) ,Y( 501 , SI 5 0) .A ( 50) , B( 5>0) .C (50) .F (50 ) ,W( 5 J) , SB( 50 ) , 
IG< 50 ). GM ( SO).?. < 50) . YI iM T< 50) ,0YDX( SO) ,02 YDXI50 ) .CiJRV ( 53 ) » R AI0( 53 ) » 

1 SOW (50) 

1 R t AO ( 5, lOoO) 

WRITE! fa. 1050) 

WRITE! fa, lOfaO) 

R EAO ( 5. 1010) N.MAX 
R tAO ( 5 . 10 20 ) ( X ( l) , l = 1 , N ) 

R EAO ( 5. I U 2 0 ) ( Y ( I ) .1=1 ,N) 

R FAD ( 5. 1020 )( Z( I ) .1 = 1. MAX) 

no io i = ? . m 

10 S ( I ) =X ( 1 )- X ( I - 1 ) 

Nll=N - ) 

IF(?.GT.U(j) GO rn 2 5 
no 20 I =2. NO 

A ( I ) = S ( I ) / 6 . 0 
B ( 1 ) =( 5 ( I ) +S( 1 + 1) ) /3. 0 

r.( i )=s( i+i >/6.o 

20 F( I ) = (Y( I + 1 ) - Y ( 1 ))/S(I+i)-!Y(l ) - Y ( 1 - 1 ) )/S(I ) 

2 5 A <N )=- .5 
B ( 1 ) = 1 . 0 
B ( N ) = 1 .0 

r.( ) ) =- . 5 

E ( I 1=0.0 
FIN) =0 . 0 
wi i )=a< i ) 

SB! 1 )=!’.( 1 ) / W ( 1 ) 

G ( 1 ) = .) . 0 

no 50 I = 2 . N 

W ( I )=5( 1 ) - A ( 1 )#SB( I - 1 ) 

S B ( I )=C< I )/W( I ) 

5 0 G< I )=( E! l)-A(l l*ii(l-ll)/WII) 
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E M ( N ) = G ( <N ) 
on 40 I = 2 • N 

K =N+ l- 1 

40 FM(K)=G(id-SH(<)*FM(K + l) 

S UM ( L) = 0. 

DO 45 1=2. N 

45 S UM I I) = S UM ( 1-1 ) + S(I M : ( Y ( I)+Y(l-1) ) / 2 . - S ( I )**3#(EM(II+EM( I — 1 ) ) / 2 4 . 
WHITE: (6.1000) N . ( X ( I ) , Y < I ) , SUM ( I ) , I =1 ,N ) 

IF(MAX.LT.l) GO TO l 
DO 90 1=1. MAX 
K =2 

I F < / ( I ) - X ( 1 ) ) 60. 5 0. 70 

5 0 Y IMT( I ) = Y( U 
GOTO 56 

60 IHI Z( l).LT.(l. l*X(l)-. l«X(2))> WRITE (S.IOOO)Z(I) 

GO TO 85 

65 I F{ Z( I ).GT ml 1. l*X( \l ) - • 1* X ( N-l ) ) ) WHITE <6,1000 )Z(I) 

K=N 

GO. TO 85 

70 IF( Z( I >— X ( X ) ) 85. /5 . 80 

75 Y INT < I > =Y ( K I 

GOTO tin 
8 0 K =K + 1 

IF(K-N) 70. 70, 65 

85 YINT(l) = EM(<-l)MX(K)-Z(I))< c *3/6./S(K.)+EM(K)*(Z(I)-X(K-l))< r *3/6. 
l/S(K)«-(Y(K)/S(X)-EM(K)*S(K)/6. )M Z ( I ) -X ( K- L ) ) + ( Y ( K-l ) / S ( < ) -EM( < - 1 ) 
2 <= S < K )/6. )* ( X(< )-Z ( I ) ) 

86 0 Y 0 X ( I ) =-EM(K- l)*( X(K)-Z 1 1 ) ) **2 /2 . 0/ S ( K ) +E M IK) * ( X ( K- 1 ) -Z ( I ) ) **2/ 2. 
10/S(K ) +< Y< K )-Y(.<- l ) ) / S(K) -(EM( K.) -EM( K-l) )*S(K)/5.0 

D2Y OX ( I ) = EMI K- I )* (X( K)-Z (I ) ) /S(K.)*EM(K) * ( Z ( I) -X( K-l) ) ZS( K ) 
CURV(l) = i)2Yi)X(I)/( l.+OYOX(I ) **2)**1.5 
R AD( I ) = 1 . /(', UK V ( 1 ) 

90 CONTINUE 

WRITE ( 6, 10 40) MAX , I Z ( I ) » YI NT( I ) » 0 YOX ( I ) .02Y0X ( I) , CURV ( 1) , R AD( l > , 

1 I = l . M A X ) 

GO TO 1 

1000 FORMAT ( 17 HI OUT IF RANGE X =GL4.6) 

1010 FORMAT ( 101 5) 

10 20 FORMAT ( 8GI0.8 ) 

1030 FORMAT ( 49HL LIST OF SPLINE COORDINATES (INPUT) (N = .12, 

1 1H) , 10X, 12 H VALUE l')F TH = / 1 5 X . I 8 X . I 9 X , 1 H Y , 2 5 X . I 5 H I NTE GR AL OF Y*D< 

2 /< 5X. 2G20.8.G32. 8) ) 

1040 FORMAT ( 1H1.5X.76HCOOROINATES. DERIVATIVES AND CURVATJRES AT SELEC 
IT FD PUNTS (NO. OF POINTS = . I 2 , 1H ) / 1 2 X , 88 X ( I NPUT ) , 8 X , 1 3 HY ( I NT E RP3 
2L AT ED) . 10X, 5HQY/0X ,14X . 7HD2Y/J X2 . 1 2 X . 9HC UR V A T URE ,6X, l 8 HR AD I US OF C 
3URVATU8E/ ( 5X.6G20. 8) ) 

1050 FORMAT ( 1H l ) 

1060 FORMAT (800 
1 


) 



FORTRAN Dictionary 


A Array of below diagonal terms of coefficient matrix 

B Array of diagonal terms of coefficient matrix 

C Array of diagonal terms of coefficient matrix 

CURV Array of curvatures at interpolation points 

DYDX Array of derivatives at interpolation points 

D2YDX Array of second derivatives at interpolation points 

EM Array of second derivatives at spline points 

F Array of constants in matrix equation . 

G Temporary array used in solving matrix equation 

I DO index 


K DO index 

MAX Number of interpolation points 

N Number of spline points 

NO N - 1 

RAD Array of radii of curvature at spline points 

S Array of length of intervals between spline points 

SB Temporary array used in solving matrix equation 

W Temporary array used in solving matrix equation 

X Array of x -coordinates of spline points 

Y Array of y -coordinates of spline points 

YINT Array of y-coordinates of interpolation points 

Z Array of x-coordinates of interpolation points 


Lewis Research Center, 

National Aeronautics and Space Administration, 
Cleveland, Ohio, September 25, 1968, 
126-15-02-31-22. 
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